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ABSTRACT 

We present a numerical scheme for modelling unresolved turbulence in cosmological adaptive mesh 
refinement codes. As a first application, we study the evolution of turbulence in the intra-cluster 
medium and in the core of a galaxy cluster. Simulations with and without subgrid scale model are 
compared in detail. Since the flow in the ICM is subsonic, the global turbulent energy contribution 
at the unresolved length scales is smaller than 1% of the internal energy. We find that the production 
of turbulence is closely correlated with merger events occurring in the cluster environment, and its 
dissipation locally affects the cluster energy budget. Because of this additional source of dissipation, 
the core temperature is larger and the density is smaller in the presence of subgrid scale turbulence 
than in the standard adiabatic run, resulting in a higher entropy core value. 

Subject headings: galaxies: clusters: general — hydrodynamics — methods: numerical — turbulence 
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1. INTRODUCTION 

Simulations of cosmological structure formation often 
share two important attributes. First, the ubiquitous 
presence of spatially localized features such as shocks, 
clumps, or composition discontinuities that need to be 
numerically resolved or at least adequately modeled; 
and second, moderate or large Reynolds numbers of the 
baryonic component indicating that fully developed, i.e. 
space-filling turbulence is responsible for the mixing and 
dissipation properties of the gas. Despite great advances 
in computational fluid dynamics, an accurate handling 
of both aspects has so far proven to be very difficult, 
because dedicated numerical techniques seem to be mu- 
tually incompatible. 

The most powerful technique for grid-based solvers to 
resolve localized and anisotropic structures in a flow is 
adaptive mesh refineme nt (AMR) (|Berger fc 01igej|1984l : 
iBerger fc Colellal Il989l ). This technique has proven to 
be very well s uited for several astrophysical problems 

orman|[2005l ). However, in the case of astrophysically 
relevant Reynolds numbers even with AMR we cannot 
resolve all the relevant length s cales down to the dissi- 
pative one (jSchmidt et al.ll2006D . Even if this condition 
may be achieved in regions of maximum refinement, as 
is possibly the case in the core of galaxy clusters (where 
the effective viscosity is still the subject of debate), tur- 
bulence from coarser areas of the grid continuously flows 
into these regions without being properly accounted for. 

In engineering applications as well as other fields 
of computational fluid dynamics, subgrid scale (SGS) 
models have been developed in order to mimic the 
influence of unresolved turbulence on the resolved 
scales. This technique is often referr ed to as Large 
Eddy Simulations (LES) (Lcsicur fc Meteii Il996l) . In 
astrophysics, SGS models have already been exten- 
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sively us ed in simulati ons of Ty pe la superno va ex- 
plosions (iNiemever fc H iUebrand fl 119951: Ul eineck e et al.l 
|20Q^ |R6Eke_^Hillebrandt 2005: 'Ro pke et al.ll2067b . In 
this framework, Schmidt ct al. (200^ presented a for- 
mulation o f SGS m o dels b ased on the filtering ap- 
proach of iGermanol ()1992l ). Other applications of 
SGS mod els in astro physical problems have been pro- 
posed by iPope et all (p008l) and, in an approach spe- 



cially designed for Rayle i gh- Ta ylor-driven turbulence, by 
IScannapieco fc Brtiggenl (j200 8). 

In this paper, a numerical method that combines LES 
and AMR for the study of astrophysical turbulent flows 
will be presented. We will refer to this new tool as fear- 
less (Fluid mEchanics with Adaptively Refined Large 
Eddy Simulations). With the combined use of grid re- 
finement and SGS model, fearless is very suitable for 
simulations of intermittent turbulent flows in clumped 
media. 

The formation and evolution of the cosmological large 
scale structure is a typical case of turbulence genera- 
tion in a strongly clumped medium. The concordance 
model of cosmological structure formation explains the 
formation of clusters through a hier archical sequen ce of 
mergers of lower-mass systems (e.g. lOstriked 119931 ) . In 
particular, mergers of subhalos play a fundamental role 
in determining the structure and dynamics of massive 
clusters of galaxies. Furthermore, it is known that major 
mergers induce temperature inhomogeneities and bulk 
motions with velocities of the o rder of 1000 km s~ ^ in th e 
intra-cluster medium (ICM) ("Norm an fc BrvanI Il999( l. 
This results in complex hydrodynamic flows where most 
of the kinetic energy is quickly dissipated to heat by 
shocks, but some part may in principle also excite long- 
lasting turbulent gas motions. Beside s merger processes, 
it is also known that galactic motions (Brcgman fc Davi 
1989: Kim 20 03) and A GN outflows (Heinz et al. 200 
Siiacki fc Springe]|l2006D can stir the ICM. 



The problem of the turbulent state of the ICM is 
still controversial, both from the theoretical point of 
view of co nstraining the kinematic viscosity of th e fluid 
(iRevnold s et al. 2005; Narayan fc Medvedev 2001; Jones! 
|2008|). and from the observational side, since a direct ob- 
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servation of turbulent emission-line broad ening is beyond 
the r e ach of current X-ray observ a tories (ISunvaey et al.l 
20031 llnogamov fc Sunwi^ 120031: iBriiggen et all 120051 
Dolag et al.n 2005: Rcbusco eTail l2008f) Nonetheless, 



some indirect ways of investigati ng turbulence in clus- 
ters g ave encouraging results (see llapichino fc Niemeyeil 
|2008| for an overview) and call for a better theoretical 
understanding of the problem. A clear example for the 
relevance of cluster turbulence for precision cosmology is 
provided by recent results that demonstrate the sensitiv- 
ity of hydrostatic mass estimates on as sumptions about 
the level of turbulence ()Lau et al.ll2009| ). 

In nu merical simulation s of merging cl usters 
dSchindlc r fc Muelleil 119931: iRoettiger et alj Il997t 
Rickey fc Saraz in' '2001' 'Fui ita et al.l l2004allbl : iTakizawal 
2005at llapichino et al. 2008[), it has been shown that 
infalling subclusters generate a laminar bulk flow 
but inject turbulent motions via Kelvin-Helmholtz 
instabilities at the interfaces between the bulk flows 
and the primary cluster gas. Such eddies redistribute 
the energy of the merger through the cluster volume 
and decay into turbulent velocity fields, eventually 
developing a turbulent cascade with a spectrum of 
fluctuation s expected to be close to a Kolmogorov 
spectrum (jPolag et al.l I2005D . Numerical simulations 
focused on the role of turbulence in astrophysical flows in 
general, and especially for clusters, have been restricted 
to measuring passively statistical qua ntities like velocity 
dispersion from simulat ion data (e.g. iNorman fc Brvaiil 
119991: iDolag et al1l2005f ). The active role of small scale 
velocity fluctuations on the large scale flow has not been 
taken into account so far. 

A previous attempt of modelling turbulence in hydro- 
dynamical simulation of cluster formatio n has been per- 
formed by llapichino fc Niemeverl (|2008f ). In that work, 
the authors focused on better deflnitions of the AMR 
criteria for reflning the computational g rid where and 
when the flow in the ICM was turbulent (| Schmidt et al.l 
l2009Hlapichino et al.|[2008l ). Though useful, this numeri- 
cal strategy can follow only a narrow range of large length 
scales along the turbulent cascade, being the Kolmogorov 
length scale for turbulent dissipation much lower than the 
spatial resolution. Besides this theoretical shortcoming, 
also numerically it is questionable whether the mixing 
forced at the mesh le ngth scale corre ctly represents the 
physics of turbulence (jMitchell et all 2009). 

These arguments motivate the application of fearless 
to cluster simulations as a more consistent approach. We 
show below that the additional degree of freedom given 
by the local turbulence intensity on unresolved scales has 
a measurable impact on the features of the ICM. In ad- 
dition to the direct dynamical coupling to the resolved 
fluid equations, the ability to separate unresolved kinetic 
energy from thermal energy allows a more accurate com- 
putation of the local temperature and entropy than with- 
out the subgrid scale model. 

This work is structured as follows: in ^ the formal- 
ism of the subgrid scale model and of fearless is intro- 
duced. Some numerical tests and consistency checks are 
presented in Sj3l and the setup of the galaxy cluster sim- 
ulations is described in f|4l The results are presented in 
^and discussed in ^\ where our conclusions are drawn. 



2.1. Germano decomposition 

The dynamics of a compressible, viscous, self- 
gravitating fluid with with density p{ri,t), momentum 
density pvi{ri, t) and total energy density pe{ri, t) at spa- 
tial position (ri,r2,r3) is given by the following set of 
equations: 



d d , , 
^P+^(-.p) = 0, 
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dri dr-j 



- P9i: 



(1) 
(2) 



where p is the pressure, gi the gravitational acceleration 
and cr'. the viscous stress tensor. Note that the Einstein 
sum conventio n applies to repeated indices. 

As shown bv l Schmidt et al.l ()2006[ ) , these equations can 
be decomposed into large-scale (resolved) and small-scale 
(unresolved) parts using the filter formalism proposed by 
iGermanol fl992. ) in terms of density- weighted quantities'*. 
By means of filtering, any field quantity a can be split 
into a smoothed part (a) and a fluctuating part a', where 
(a) varies only at scales greater than the prescribed filter 
length. We define densi ty weighted filtered quantities 
according to iFavrd (|1969D by 



(pa) = (p)a => a - 



(4) 



Following ISchmidt et al.l (|2006D . filtered equations for 
compressible fluid dynamics can be derived: 
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— (P> + T5-iij(P> = 0, 
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+ {p){X + e)-v,—i-{v„Vj) (7) 
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ar-j 

where we introduced the total resolved energy Cres = 
Gint + \'ViVi, being eint the filtered internal energy, and 
the generalized moments which are generically defined 

by 



f (a, b) = {pab) — {p)db 
T{a, b, c) = (pabc) — {p)abc 

— drib, c) — 6f (a, c) — cf (a, 6) 
f (a, 6, c, d) = ... 



(8) 

(9) 
(10) 



for Favre-filtered quantities a, 6, c etc. Germano inter- 
preted the trace of f{vi,Vj){p) as the squared velocity 
fluctuation, q'^ :— T{vi,Vi)/ (p). The evolution of the cor- 

■* For a review, see lRopke fc Schmidt) ||2009|1 . 
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responding turbulent energy, et 



, is given by 



dr 



-Vj{p)et 



) + S + r- (p>(A + e) 



where 
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(9r, 



T{Vj,Vi,Vi) - ^J, + I 



r : 

(P>A: 



-Tijdvi/drj 
+r{vi,gi) 
9 . 
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dri 



dri 



d 



d_ 



and 



-p = {Vip) - Vi{p) 

-K = {Via'ij) - Viia'ij) 



(11) 



(12) 

(13) 
(14) 

(15) 
(16) 



(17) 
(18) 



The explicit forms of the quantities D, A, e, F and f (w^, vj) 
are unknown and have to be modeled in terms of closure 
relations, i. e., functions of the filtered flow quantities 
(or their derivatives) and the turbulent energy et. The 
closures for all these terms represent the SGS model. 

2.2. Subgrid scale closures 

In the following we consider a simplified set of equa- 
tions to model the influence of the turbulent small scale 
(SGS) motions on the numerically resolved scales £> £a, 
neglecting the influence of the viscous stress tensor (cr^^ ) 
(which is a very good approximation for high Reynolds 
numbers) and the turbulent transport of heat given by 
the divergence of T('!;j,eint) in equations ([6]) and ([7]). 
Moreover, gravitational effects on unresolved scales are 
neglected, i. e., we set F = in equation (fTT|) for the 
turbulent energy. For the terms (fT^ . and we 
adopt SGS closures that have been applied in large eddy 
simulatio ns of incompressible t urbulence. The numerical 
study bv lSchmidt et al.l (|2006D demonstrated that these 
closures can be carried over to transonic turbulence, for 
which the unresolved turbulent velocity fluctuations are 
small compared to the speed of sound. Additionally , 
we utilize the pressure-dilatation model of lSarkaj ()1992f ) 
in order to account for moderate compressibility effects. 
Since we concentrate on the dynamics of the gas in the 
ICM, where the Mach numbers may locally approach 
unity, the SGS closures outlined subsequently serve as 
a reasonable approximation. In supersonic flow regions, 
on the other hand, the SGS model is deactivated in order 
to maintain stability (see ^12.4^ . 

The flux of kinetic energy from resolved scales toward 
subgrid scales, i. e., the rate of turbulent energy produc- 
tion, is given by the contraction of the turbulent stress 
tensor and the Jacobian of the resolved velocity field. 
Since Tkk = (p)?^; we split the tensor in a symmetric 
trace- free part T*j and a diagonal part: 



(19) 



The mod el for f is based on the turbulent viscosity 
hypothesis (jBoussinesd 118771) . which means that f,* is 
assumed to be of the same form as the stress tensor a',- 



of a Newtonian fluid. Hence, 

= —2r]tSij 

with a turbulent dynamic viscosity rjt 
{pjCiylAq and 



(20) 



dri 



3 ' drk 



(21) 



The turbulence production term is therefore modeled as 

1, 
3' 



(22) 



We set = 0.05 (lSagautJl200d ). 

The SGS transport of turbulent energy feauation[T2|) is 
modeled by a gradient-diffusion hypothesis, stating that 
the non-linear ter m is proportio nal to the turbulent ve- 
locity gradient (|Sagautil2006( ) 



d 8 
dri dri 



(23) 
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The diffusion coefficient ha s been calibrated to Ci 
by numerical experiments (jSchmidt et al.|[2006r ). 

For sufficiently high Reynolds numbers, viscous energy 
dissipation (equation [T6|) becomes entirely an SGS effect. 
The most simple expression that can be built from the 
characteristic turbulent velocity and length scale for dis- 
sipation is 



(a 



(24) 



For our simulations we set = 0.5 ()Sagautil2006r i. 

The effect of unresolved pressure fluctuations in com- 
pressible turbulence is described by the A term (equa- 
tion [151). A simple closure for subsonic turbulent flow is 
(Deardora.1973.1 



A = C^q^ 



(25) 



where Ca = (iFurebv et al.lllQQl . iSarkad (fTOOl per- 
formed simulations of simple compressible flows and in- 
vestigated the influence of the mean Mach number of 
the flow on the turbulent dissipation e and the pressure 
dilatation A. Based on this analysis he suggested dif- 
ferent models for these terms, which we will describe in 
the following sections. These modifications have been 
proven to yield good result s for transonic turbulence 
()Shvv fc K rishnamurtv 1997^. 

As a major effect of_compressibility from direct numer- 
ical simulation, ISarkarl (|i99^ identified that the growth 
rate of kinetic energy decreases if the initial turbulent 
Mach number increases. This means that the dissipation 
of kinetic energy (and, therefore, of the turbulent energy) 
increases with the turbulent M ach number M t = q/cs, 
where Cg is the speed of sound. ' SarkaH ()1992D suggested 
to account for this effect by using 



e = C,|-(l + aiM2) 
'A 



(26) 



with ai = 0.5 as a model for the dissipation of turbulent 
energy. 

Based on a decomposition of all variables of the equa- 
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tion for instantaneous pressure 
92 



dridr 



■{pViVj 



(27) 



into a mean and a fluctuating part and comparisons with 
direct numer ical simulations of simple compressible flows 
iSarkarl ()1992f ) proposed a different model for the pressure 
dilatation 



A = 02Mtf * 



'*7 



'a 



■ 8a4M^pt 



dru 



(28) 



with a2 — 0.15,a3 — 0.2 obtained from a curve fit of 
the model with direct num erical simulation (DNS). Un- 
fortunately, iSarkail (|1992f ) does not specify a value for 
a4, so there i s some confusion in the l iterat ure about it. 
For example, IShvv &: Krishnamurtvl (|1997f) set 04 = 
and still found the Sarkar model in good agreement with 
their DNS simulation. In this work, we adopt ot^ = a\/2. 
With this choice, the effective production of turbulent en- 
ergy vanishes for a turbulent Mach number Mt = l/a2- 
In the following we will sometimes refer to the "Sarkar 
SGS" when the equations (|^ and (|^ are used in the 
model. 



2.3. Filtered equations in comoving coordinates 

Simulations with a comoving cosmological background 
require a formulation of the filtered fiuid dynamical equa- 
tions in comoving coordinates. Applying the Germano 
decomposition f H2.ip in a comoving coordinate system 
with spatially homogeneous scale factor a(i), we obtain^ 



dt 



a dxj 



{p}Ui + --^Uj{p}Ui 



at 



a dxi 



l_d_ 

a dxi 
l_d_ 

a dx. 



(p) + {p)9i 



(29) 



(30) 



■T{Ui,Uj) {p)Ui 



d_ 

at 



■^(p)eres -I- --^Uj{p)ereB = - - -^Ui{p) - - {p)u,g* 

a axj . j_ 



1 d 

a dxi 



1 



-{{P)e 



res -r -{p)UiUi + (p>) 



+ (p)(A + e) - -Ui^—T{ui,Uj) 



dxi 



(31) 



^ (p>et + -^Uj{p)et = D + r- (p>(A + e) 



dt 



a dxj 



T{uj,Ui)- — Ui - 2-{p)et 

a ax-i a 



(32) 



W ith respect to the non-comoving equations listed in 
'> \2.1[ the only term to implement additionally is the last 
one on the right-hand side of equation ([32|l . Further- 
more, SGS closures have to be expressed in terms of the 
Jacobian of the velocity in comoving coordinates. 



Jij — — Vj — — — Uj H Oi- 
avi a oxi a 



(33) 



In particular, the trace- free rate of strain tensor in co- 

^ The comoving density p = a^p and the comoving pressure 
p = a^p are introduced to shorten the equations. 



moving coordinates is given by 



(34) 



2.4. Limits of the SGS model 



Numerical difficulties result from constraints on the va- 
lidity of SGS closures. In particular, the turbulent vis- 
cosity hypothesis expressed by equation ([22]) was devised 
to account for the production of turbulence by shear in 
moderately compressible fiow. This is typically encoun- 
tered in the dense, central regions of galaxy clusters. 
However, the surrounding low-density gas can be acceler- 
ated very quickly in the gravitational field of the cluster. 
Moreover, high velocity gradients are encountered in the 
vicinity of shocks which are produced by gas accretion 
onto filaments, sheets, and halos as well as by the merg- 
ing of substructures. In order to inhibit unphysical pro- 
duction of turbulent energy by these mechanisms, which 
are not accommodated in the present formulation of the 
SGS model, several numerical safeguarding mechanisms 

have been introduced. 

First of all, we implemented a simple shock detector 
which identifies strong negative divergence. A cell is 
marked if the velocity jump corresponding to the neg- 
ative divergence becomes greater than the sound speed 
across the cell width Ia: 



dri I A 



(35) 



In cells satisfying the above criterion, the source and 
transport terms of the SGS model (equations and [^5)1 
are disabled. The turbulent energy is only advected in 
these cells, and no coupling to the the velocity and the 
resolved energy is applied. 

Besides the previous check, an additional constraint 
is imposed on the magnitude of the turbulent energy, 
via the turbulent Mach number. Basically, the SGS 
model breaks down once Mt becomes large compared 
to unity, therefore a threshold for this quantity is set 
to Mt^inax = V^. This value for the maximal turbulent 
Mach number is motivated by the theory of isothermal 
turbulence, where the effective gas pressure can be ex- 
pressed as 



Poff = PCs 



-p(j2 = 7<,ff pc^ 



(36) 



and, consequently, Pcff is limited to the adiabatic value 
|pCg. We verified that this threshold does not harm 
our cluster simulations, because in the hotter gas phases 
(T > 10^ K) turbulence is largely subsonic, and the 
threshold is rarely reached. 

A supplementary low temperature cutoff ensures that 
the sound speed does not drop to excessively low values, 
which occur in cosmological simulations especially in the 
low-density voids. We set the lower limit of the temper- 
ature to Tmin = 10 K. This threshold ensures numerical 
stability and does not affect the baryon physics apprecia- 
bly, apart from possibly making the shocks on accreting 
structures weaker. 

2.5. Combining AMR and LES 

In Large Eddy Simulations (LES), the filtered equa- 
tions (|29ll32p are solved using an SGS model as out- 
lined in the previous Section. However, the closure rela- 



Adaptively refined large eddy simulations of clusters 



5 



tions we use and, in fact, the very concept of SGS tur- 
bulence energy only applies if the velocity fluctuations 
on subgrid scales are nearly isotropic. This limits the 
LES methodology to flows where all anisotropics stem- 
ming from large scale features, like boundary conditions 
or external forces, can be resolved. In the fearless 
method, the grid resolution is locally adjusted by 
adaptive mesh refinement (AMR) in order to ensure that 
the anisotropic, energy-containing scales are resolved ev- 
erywhere. On the other hand, it is assumed that turbu- 
lence is asymptotically isotropic on length scales compa- 
rable to or less than ^a- It is very difficult to justify the 
latter assumption a priori, because there are no refine- 
ment criteria that would guarantee asymptotic isotropy 
on the smallest resolved length scales. By careful analysis 
of simulation results, however, one can gain confidence 
that AMR resolves turbulent regions appropriately. 

As an infrastructure for the implemen tation of FEAR- 
LESS, we chose Enzo v. 1.0 (jO'Shea et al.|[20 05l. an AMR, 
grid-based hybrid (hydrodynamics plus N-Body) code 
based on the PPM solver (Colella & Woodward 198|) 
and especially designed for simulations of cosmological 
structure formation. When a grid location is flagged 
for refinement in Enzo, a new finer grid is created, and 
the cell values on the finer grid are generated by in- 
terpolating them from the coarser grid using a conser- 
vative interpolation scheme. At each timestep of the 
coarse grid, the values from the fine grid are averaged 
and the values computed on the coarse grid (in the re- 
gion where fine and coarse grid overlap) are replaced. 
However, this approach does not account for the inher- 
ent scale-dependenc e of the turbulent energy. Assu ming 
Kolmogorov scaling (|Kolmogorovlll94ll lFrischllT995D , the 
turbulent energies at two different levels of refinement 
with cell size and Ia,2, respectively, are statistically 
related by 

'^^i^C^yr (37) 

et,2 (?| V'A,2/ 

Using this scaling relation, we implemented a simple 
algorithm to adjust the turbulent energy budget when 
grids are refined or derefined. The following procedure 
is used once a grid is refined: 

1. Interpolate the values from the coarse to the fine 
grid using the standard interpolation scheme from 
Enzo. 

2. On the finer grid, correct the values of the ve- 
locity components, i'i, and the turbulent energy, 
et = ^q^, as follows 

= ^' Jl + ^fl-rX'/'), (38) 

V ekin ^ ^ 

e't = etr--^/^ , (39) 

where ekin is the resolved kinetic energy, ta is 
the refinement factor of the mesh, and the primed 
quantities are the final values on the fine grid. The 
resolved energy is adjusted such that the sum of 
resolved energy and turbulent energy remains con- 
served. 

Apart from adjusting the energy budget, the resolved 
flow should feature velocity fluctuations on length scales 
smaller than the cutoff length of the parent grid if a 



refined grid is generated. To address this problem we 
observe that the smallest pre-existing eddies that are in- 
herited from the parent grid will produce new eddies of 
smaller size within a turn-over time. Although this im- 
plies a small delay because of the higher time resolution 
of the refined grid, the flow will rapidly adjust itself to 
the new grid resolution. 
For grid derefinement we reverse this procedure: 

1. Average the values from the fine grid and replace 
the corresponding values on the coarse grid 

2. In the regions of the coarse grid covered by flner 
grids, correct the velocity components and the tur- 
bulent energy: 

e[ = etri/^ (40) 

<-''f^IW^)- '''' 

Here, primed quantities denote the flnal values on 
the coarse grid. The resolved energy is adjusted to 
maintain energy conservation and a positive kinetic 
energy. 

3. NUMERICAL TESTS 

We applied two consistency tests of the SGS model in 
simulations of forced isotropic turbulence in a periodic 
box. First, energy conservation was checked in adiabatic 
turbulence simulations and, second, the scaling of the 
turbulent energy over several levels of resolution was in- 
vestigated for isothermal turbulence. To simulate driven 
turbulent flow, a random forcing mechani sm ba sed on the 
Ornstein-Uhlenbeck process was applied (jSchmidt 2004). 
This process generates a solenoidal (i. e., divergence-free) 
stochastic force field which accelerates the fluid at large 
length scales I ^ Iq, where Iq is the size of the computa- 
tional box. The strength of the force field is character- 
ized by a forcing Mach number Mf . The Mach numbers 
of the fiow becomes comparable to Mf once the forcing 
has been applied over a period of time that is defined by 
the integral time tint — lo/MfCg. 

3.1. Energy conservation 

For global energy conservation it turned out to be im- 
portant to compute the turbulent stress term in equation 
([6|) indirectly as 

d d d 

ViT;—rivi,Vj) = ——ViT{vi,Vj)-T(vj,Vi)^—Vi , (42) 

The reason behind it is that only by using this rearrange- 
ment of the terms we can ensure that we do not introduce 
small numerical errors which would violate the local sum 

— ViT{Vi,Vj) = Vi—T(Vi,Vj) + T{Vj,Vi)—Vi. (43) 
OTj OTj OTj 

leading to a big error in global energy conservation. 

As a testing case, we run a LES of driven turbulence 
as outlined above on a static grid of 256'^ grid points 
and periodic boundary conditions. The adiabatic index 
7 = f and Mf = 0.68. 

Figure [T] shows the typical time development of the 
mass- weighted mean energies in our simulation including 
the energy Cf injected into the system by random forcing. 
It is evident from the curve of the turbulent energy that 
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256 , Ma = 0.7, Basic model. 



256- . Ma = 0.7 
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Fig. 1. — Time evolution of mass-weighted energy averages in 
the driven turbulence simulation. The different energy components 
are indicated by colors. 
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Fig. 3. — Time evolution of the energy differences (red line: 
kinetic energy; green line: internal energy) between simulations 
with and without the SGS model (basic version), compared to the 
evolution of the turbulent energy (blue line). 
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Fig. 2. — Relative error of the total energy in the driven turbu- 
lence simulation. The different lines indicate simulations run with 
different versions of the SGS model, or without it, as shown in the 
legend. 

after one integral time scale, our simulation reaches an 
equilibrium between production and dissipation of tur- 
bulent energy. 
In Fig. [21 we plot the time development of the relative 

error of the mean total energy, defined as the sum 

of the mass-weighted means of internal energy, kinetic 
energy, turbulent energy minus the injected energy by 
the forcing, 



eint + fikin -I- et - ef 



(44) 



where e ~ ^r^- It demonstrates that with our basic 
model, the relative error in energy is comparable to the 
error without SGS model, and is around 1%. The energy 
conservation of the model using the Sarkar modifications 
is equally good. 

It is also instructive to plot the difference between the 
energy contributions (internal and kinetic energy) in the 
simulations with and without the SGS model. These dif- 
ferences are shown in Fig. [S] One can conclude from this 
figure that, at the beginning of the turbulent driving, the 
turbulent energy produced in our simulation with SGS 
model is found in the kinetic energy of the simulation 
without SGS model, in contrast, from t = 1.2 tint on, 
most of the turbulent energy can be found in the internal 
energy of the simulation without SGS model. Turbulent 
energy can therefore be interpreted as a kind of buffer 



AMR level 
AMR level 1 
AMR level 2 
AMR level 3 




Fig. 4. — Thick lines: mean mass-weighted turbulent energy for 
each level of the AMR simulation, using our procedure of transfer- 
ring turbulent energy at grid refinement/derefinement. Thin lines: 
the corresponding evolution of turbulent energy of the static grid 
simulations. The colors indicate the AMR level or the static grid 
resolution. 

which prevents the kinetic energy in our simulation to 
be converted instantly into thermal energy. 

3.2. Scaling of turbulent energy 

A necessary condition for the validity of the turbu- 
lent energy transfer algorithms explained in !j2.5l is that 
an AMR simulation should approximately reproduce the 
results of static grid simulations corresponding to the dif- 
ferent levels of refinement. To test this, we compared an 
AMR simulation of driven turbulence with a 32"^ root 
grid resolution and three additional levels (with a refine- 
ment factor of 2 between each level) to three static grid 
simulations with resolutions of 32'^, 64'^, 128"^ and 256'^. 
In order to allow for the comparison of averaged quanti- 
ties, refinement of the entire domain was enforced at all 
levels of the AMR simulation. In order to reach a statis- 
tically stationary root-mean-square (rms) Mach number, 
we ran the simulations for nearly isothermal gas with 
7 = 1.01 (for adiabatic turbulence, after an initial rise, 
the rms Mach number gradually decreases with time be- 
cause of the dissipative heating of the gas). We used 
a supersonic forcing Mach number, Mf = 2.7, to check 
whether a consistent turbulent energy budget could be 
achieved for highly compressible turbulence, albeit the 
scaling relations ([37]) of incompressible turbulence were 
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Fig. 5. — Thick lines: mean mass-weighted turbulent energy 
for each level of the AMR simulation without transferring tur- 
bulent energy at grid refinement/derefinement. Thin lines: the 
corresponding evolution of turbulent energy of the static grid sim- 
ulations. The colors indicate the AMR level or the static grid 
resolution. 

utilized. 

The results of this consistency check can be seen in 
Fig. m We observe that the time development of the 
mean turbulent energy is very similar on the different 
levels of the AMR simulation compared to the static grid 
simulations, except for some deviations at the root level. 
On the other hand, comparing these results to a simula- 
tion without correcting turbulent energy at grid refine- 
ment / derefinement (Fig. [5|) , it is evident that the scaling 
of the turbulent energy in the latter case is inconsistent. 

4. DETAILS OF THE CLUSTER SIMULATIONS 
4.1. Simulation setup 

We performe d simulations of cluste r forni ation with 
Enzo, following 'lapichino fc Niemeverl (|2008l). We will 
compare two runs: one of them was done with the pub- 
lic version of Enzo (without SGS model) , and the other 
with FEARLESS implemented, in its version including the 
Sarkar correction (equations [25] and The simula- 

tions were done using a flat ACDM background cosmol- 
ogy with a dark energy density — 0.7, a total (includ- 
ing baryonic and dark matter) matter density flm = 0.3, 
a baryonic matter density = 0.04, the Hubble pa- 
rameter set to h = 0.7, the mass fluctuation amplitude 
as = 0.9, and the scalar spectral index n = 1. Both 
simulations were started with th e same initial con ditions 
at redshift Zjni = 60, using the lEisenstein fc Hul (|l999f ) 
transfer function, and evolved to z = 0. The simulations 
are adiabatic with a heat capacity ratio 7 = 5/3 assum- 
ing a fully ionized gas with a mean molecular weight 
= 0.6 u. Cooling physics, magnetic fields, feedback, 
and transport processes are neglected. 

The simulation box has a comoving size of 
128 Mpc h~^. It is resolved with a root grid (level I = 0) 
of 128^ cells and 128^ N-body particles. A static child 
grid (Z = 1) is nested inside the root grid with a size 
of 64 Mpc h-\ 128^ cells and 128^ N-body particles. 
The mass of each particle in this grid is 9 x 10^ M© h~^. 
Inside this grid, in a volume of 38.4 Mpc h~^, adaptive 
grid refinement from level Z = 2 to Z = 7 is enabled 
us ing an overdensity refineme nt criterion as described 
in llapichino fc Niemeyeil ()2008l ) with an overdensity fac- 
tor / = 4.0. The refinement factor between two levels 
was set to TA —2, allowing for an effective resolution of 



7.8 kpc h-\ 

The static and dynamically refined grids were nested 
around the place of formation of a g alaxy cluster, identi- 
fied using the HOP algorithm (Eise nstein fc Hutl ll9981. 
Since the realization of the initial c onditions was chosen 
identical to lapichino fc Niemeveil (|2008. ) . this study is 
based on the same cluster analyzed in that work. The 
cluster has a virial mass of Mvir = 5.95 x lO^'^M© h~^ 
and a virial radius of i?vir = 1.37 Mpc h~^. 

4.2. Local Kolmogorov scaling 

In static grid simulations one often chooses to use the 
grid resolution Za as characteristic length scale to com- 
pute a characteristic velocity or eddy turnover time for 
this scale. However, in an AMR code it is not trivial 
to compute the turbulent velocity qi associated with a 
characteristic length scale I = Ia, since I a varies in time 
and space. To circumvent this difficulty, we assume that 
below the grid resolution turbulent velocity locally scales 
according to Kolmogorov 

q{l) ~ Z^/^ . (45) 

We thereby assume that locally a Kolmogorov- like energy 
cascade sets in, at a length scale given by the resolution of 
the grid at that position. This local hypothesis holds here 
only for the analysis of our simulations, and is similar 
to the assumption done in W2.5\ for managing the grid 
refinement and dereflnemnt. 

As a characteristic scale of our analysis, we choose 
the length scale of our highest resolved regions, which 
is Imin = 7.8 kpc h~^. The turbulent velocity in the most 
finely resolved regions can be computed directly from 
the values of the turbulent energy q{l) — ^/2et on the 
grid; the turbulent velocity in less finely refined regions 
is scaled down according to our local Kolmogorov hy- 
pothesis as 

q{lmin) = q{lA)(^j^j ■ (46) 
5. RESULTS 

5.1. Turbulent energy scaling in the cluster simulation 

In ^ 13.21 we studied the temporal evolution of the tur- 
bulent energy at different resolutions in a simulation of 
driven turbulence. In this section, we repeat this analysis 
for our FEARLESS cluster simulation. Figure [H] shows the 
evolution of the mass-weighted mean turbulent energy 
for every level of our AMR simulation. We see from the 
plot that the turbulent energy on the higher AMR lev- 
els I (meaning at smaller scales) is higher at early times 
(2 Gyr < t < 6 Gyr). Later this picture changes, but not 
completely. For example, the turbulent energy at Z = 4 
stays above the turbulent energy at Z = 3 throughout the 
simulation. The magnitude of et along the AMR levels 
suggests to locate the turbulence injection length scale 
between 125 and 250 kpc h~^, corresponding to the ef- 
fective resolutions of levels 3 and 4. This is only a rough 
qualitative estimate, but nev ertheless in agreement wit h 
theoretical expectations (e.g. ISubramanian et al.ll2006h . 

Particularly noteworthy are the turbulent energy fluc- 
tuations on smaller scales at the time 2 Gyr < t < 6 Gyr, 
corresponding to a redshift z = 3—1. We can inter- 
pret these large fluctuations as evidence for violent ma- 
jor mergers that happen at that time, producing turbu- 
lent energy which is then dissipated into internal energy. 
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Fig. 6. — Time evolution of the mean turbulent energy for each 
level of refinement. This analysis has been performed on a test run 
identical to that described in i|4.1l with only six AMR levels. 

heating up the cluster gas. However, at < > 12 Gyr the 
simulation seems to globally reach some kind of stable 
state, comparable to what has been found in the driven 
turbulence simulations. 

5.2. Spatial distribution of turbulent energy 

Before performing a quantitative analysis of the cluster 
properties in the fearless run, it is useful to visually 
inspect the generation and the spatial distribution of the 
turbulent energy in the ICM and around the cluster, in 
order to compare the simulations with the theoretical 
expectations of cluster mergers. Figs. [7] and [8] present a 
time series of density and turbulent velocity slices, where 
several merger events in the cluster outskirts can be iden- 
tified. 

In the density slice at redshift z = 0.15 (Fig. [TJi), 
we can see a filament extending from the lower left to 
the upper right corner of the figure; material is falling 
onto the cluster along this structure. On both sides, 
the inflow of relatively cold gas from the filament onto 
the ICM produces a moderate increase of turbulent en- 
ergy (cf . Nagai fc Kr avtsov 2003) . From the upper right 
side there is not only a smooth inflow of matter, but 
two small clumps are approaching the cluster. During 
the simulation these two clumps merge with the main 
cluster (Figs. ^ to [8h) and one of them (on the left) is 
assimilated completely at redshift z = 0. A substructure 
approaches the cluster from the lower left corner along 
the filament (Fig. [8jz) and another one is visible only at 
z = 0.05 just to the right of the cluster core, when it 
crosses the slicing plane. 

The merging process can be followed much more easily 
in terms of production of turbulent energy, visualized 
by the turbulent velocity q = v^Set. In Fig. [7f) at z = 
0.15, a marked peak of turbulent energy in the center of 
our cluster, resulting from a former massive merger, can 
be seen. The turbulent energy produced by this merger 
declines (Figs.[7li to[8ji) and at z = it is dissipated into 
internal energy nearly completely, as confirmed by our 
further analysis in §5.41 

The two approaching clumps described above continue 
to drive turbulence in the cluster. Thereby, the left 
clump can be identified in the turbulent velocity slice at 
z = 0.15 (Fig. [7{)) as a ring-like structure, showing that 
turbulence is not produced at the center of the infalling 



TABLE 1 

Energy contributions in a sphere of r = Rvir, 

CENTERED AT THE CLUSTER CORE, AT 2 = 0. 



Quantity 


Adiabatic run 


FEARLESS run 


Etot [10^3 gj.g] 


2.6458 


2.6426 (-0.1%) 


Eint [1063 erg] 


2.1982 


2.2082 (+0.5%) 


Ekin [1062 erg] 


4.476 


4.168 (-6.9%) 


Et [10^1 erg] 




1.762 



Note. — The total energy Etot is defined as the 
sum of Eint, Ejtin and, in the FEARLESS run, Et. 
The turbulen t ene rgy reported here is not scaled as 
described in i|5.1l 



clump but at the front (behind a bow shock) and in the 
wake of the infalling material. The right clump only 
shows some turbulence production in its wake, which 
might be due to its smaller size and smaller velocity. On 
their way towards the main cluster and through its ICM, 
both clumps show a relevant production of turbulent en- 
ergy (Figs.[7Ji to[Sf)). The considerable amount of turbu- 
lent energy can even be identified after the two clumps 
have merged with the main cluster (Fig. [51i) and the 
left one is not easily visible in the corresponding density 
slice (Fig. [S];). From this point of view, the distribution 
of turbulent energy traces the local merging history of a 
galaxy cluster until it is dissipated into internal energy 
completely. 

The morphological evolution of the cluster gives a clear 
sense of the markedly local behavior of the production 
and dissipation of turbulence, which is confirmed to be 
an intermittent process in the ICM. 

5.3. Cluster energy budget 

It is extremely difficult to apply an energy analysis 
similar to that performed in H3.ll to a galaxy cluster. 
Different than a periodic box, a galaxy cluster is an open 
system, with a growth over time of negative gravitational 
potential energy. Nevertheless, a comparison of the en- 
ergy contributions of the two simulations at z = is 
useful to understand the role of the SGS model. 

The results are summarized in Table [TJ where the en- 
ergies in a sphere centered at the cluster center and 
with r — i?vir are reported. Different than elsewhere 
in this work, in Table [T] the energies are not specific, 
i.e. Etot = petot etc. This choice allows a better eval- 
uation of the energy budget but it docs not differ ap- 
preciably from the analysis of specific energies, since the 
baryon masses in the two runs agree within 1%. 

The total energy remains basically unaltered in the two 
simulations, whereas the most important change is the 
decrease of Skin in the FEARLESS run. The missing ki- 
netic energy is transferred mostly to Et, which acts as 
a buffer between the resolved kinetic energy and i?int- 
The SGS model transfers energy from Et to Ei^t either 
adiabatically (via the pressure dilatation term, equation 
[T5|) or irreversibly (via the dissipative term, eauationfTOI. 
Turbulent dissipation is thus added to the dominant nu- 
merical dissipation, resulting in a moderate increase of 
i?int , though it is less relevant than the variation of -Bkin ■ 

In both runs, the kinetic energy contribution in the 
cluster -Bkin is smaller than Ei^f The mass- weighted 
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Fig. 7. — Slices of baryon density (left-hand panels, a and c) and turbulent velocity q = \/2et scaled to 7.8 kpc h ^ (right-hand panels, b 
and d) at different redshifts z, for the FEARLESS run. The density is logarithmically color coded as overdensity with respect to the average 
baryon density in the colorbar on the left of panel a, whereas q is linearly coded in km s~^, according to the colorbar on the left of panel 
b. The overlayed contours show density. The slices show a region of 6.4 X 6.4 Mpc h~^ around the center of the main cluster followed in 
the simulation. Panels a and b refer to z = 0.15, panels c and d to 2 = 0.1. 

average of the Mach number in the cluster is about 0.6, 
in agreement with the known fact that the flow in the 
ICM is, on average, mildly subsonic. In this regime it 
is not surprising that the subgrid energy Et is about 
two orders of magnitude smaller than iSint on the global 
level. The energy contribution from unresolved scales is 
globally negligible, though locally turbulence can play a 
more significant role, as will be discussed in ^5A\ 

Finally, we note that the ratio of the turbulent produc- 
tion term S (equation I13p to the turbulent dissipation 
term e (equation I16p in the cluster core is 



TABLE 2 

Mass-weighted averages in a sphere of r = 

CENTERED AT THE CLUSTER CORE, AT 2 : 



0.07 i?vir 
= 0. 



: 0.93 



(47) 



suggesting that the turbulent flow in the ICM is glob- 
ally in a regime of near equilibrium of production and 
dissipation of turbulent energy. 

5.4. Radial profiles and local analyses 

The results of referring to global features of the 
galaxy cluster, will be complemented by a local compar- 
ison in terms of radial profiles of selected physical quan- 



Quantity 


Adiabatic run 


FEARLESS run 


S/e 




0.59 


etot [lOlf^ cm2 s-2] 


1.3781 


1.4189 


eint [10" cm^ s-2] 


1.2529 


1.3030 


ejiin [10^^ cm2 s-2] 


1.252 


1.138 


et [10" cm2 s-2] 




2.078 


Pres/(Pros + Pthcrm) [%] 


1.46 


1.52 


Drms [km S^^] 


196 


204 



titles and an analysis of the cluster core in this section. 

As a first interesting result of the comparison between 
radial profiles in the fearless simulation and in the 
standard adiabatic run, one can notice (Fig. [5]) that the 
temperature profile of the fearless run deviates slightly 
from that of the adiabatic run. This is especially ap- 
parent at the center, where T is larger in the cluster 
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Fig. 8. — Same as Fig. [3 but in this case panels a and b refer to z = 0.05, and panels c and d to z = 0. The black rectangle in panel a 
denotes the projection on the slice of a small volume, including a subclump and its wake, analysed in Table [3] 



FEARLESS 
adiabatic 



Fig. 9. — Radial profiles of mass- weighted temperature at z = 
0. The dotted line refers to the simulation without SGS model, 
whereas the solid line is for the FEARLESS simulation (in its version 
including the Sarkar corrections). 



core for r < 0.07 -Rvir, with respect to the standard run. 
Consequently (Fig. [TO|l the core in the fearless run is 
less dense, so that the ICM remains in hydrostatic equi- 
librium. The local energy budget in the cluster core is 
therefore modified by the SGS model. 




Fig. 10. — Same as Fig. [9] but showing the baryon density. 

In Table [21 we explore this feature in detail, report- 
ing the mass- weighted averages of selected variables in a 
sphere within 0.07 i?vir from the cluster center. Because 
of the adjustment of the cluster hydrostatic equilibrium, 
the mass enclosed in this sphere is significantly different 
in the two runs (it decreases by 10 % in the fearless 
run), thus it is more convenient to present specific ener- 



Adaptively refined large eddy simulations of clusters 



11 




FEARLESS 
adiabatic 



Fig. 11. — Radial profile of the turbulent contribution to the 
pressure support pt/(pt +Pthcrm)i defined in the text, for the 
FEARLESS simulation at 2 = 0. 

gies in Tabled 

First, the low value of the S/e ratio indicates that, 
at z = 0, in the cluster core region the dissipation of 
turbulence is dominant with respect to its production. 
This confirms the result of the morphological analysis 
in §5.21 that turbulence is not produced locally in the 
core by mergers a.t z < 0.15, but that it decays in this 
region. The impact of this turbulent dissipation on the 
local energy budget of the cluster core can seen from 
the comparison of the energy contributions in Table O 
Similar to the global analysis in Table [H there is a clear 
decrease of ekin, transferred both to et and eint- Both etot 
and eint are higher in the FEARLESS run, pointing to the 
existence of an energy flux from the resolved scales to the 
thermal reservoir through the turbulent buffer, leading 
to the increase of the internal energy. We interpret this 
additional energy contribution as caused by the turbulent 
dissipation introduced by the SGS model. 

In the cluster core, the energy content at the subgrid 
scales is marginal. Apparently the relative contribution 
to the total energy is even smaller than in Table [U but 
one should notice that, for consistency, in that table 
both E'kin and Et are reported according to the origi- 
nal scale separation introduced by the AMR resolution, 
and without rcscaling Et as described in %.![ In the 
cluster core the refinement level is maximum, therefore 
the unresolved part of the turbulent cascade is relatively 
smaller than elsewhere, and so is et- In Table [2] we use 
the scaled definition of et, but in the core it differs from 
the unsealed one only marginally, because in this region 
the resolution is /min almost everywhere. 

To further quantitatively appreciate the contribution 
of et to the energy budget. Fig. [11] reports the pro- 
file of the turbulent pressure support pt/ipt + Pthcrm) 
in the cluster, where the turbulent pressure is defined 
as pt — 1/3 pq^, and pthcrm is the usual thermodynami- 
cal pressure. This ratio is also equal to the ratio of the 
corresponding energies (et/(et -|- eint)). At the length 
scale of the effective spatial resolution of the simulations 
'mill = 7.8 kpc h~^, the contribution of the turbulent 
pressure (or energy) is well below 1%, although it in- 
creases at larger central distances. 

In analogy with the turbulent pressure, we define a 
"resolved pressure" 



Fig. 12. — Same as Fig. |9] but showing the mass- weighted en- 
tropy (as defined in the text). 




Fig. 13. — Radial profile of th e tu rbulent energy scaled at the 
length scale iniinj ^ described in i|4.2l for the FEARLESS simulation 
at 2 = 0. 



where the root-mean-squar e (hereafter rms) 
(jlapichino &: NiemeveH [20081 ) is defined as 



velocity 



(49) 



Here, (v) is the mass- weighted average of the velocity in 
the analysis volume. This quantity essentially probes the 
contribution of turbulent motions at length scales of the 
order of 0.07 i?vir ^ 90 kpc h~^. As shown in Table [2 
the pressure contribution at these length scales is at the 
percent level, and is slightly higher in the fearless sim- 
ulation. Interestingly, the rms velocity in the FEARLESS 
run is also somewhat larger than in the adiabatic case. 

The changes in the temperature and density profiles 
are also reflected on the entropy which is defined, as is 
customary in astrophysics, as 



K : 



(50) 



Prcs ■ 



L 2 



(48) 



with 7 = 5/3. The entropy in the cluster core is higher 
in the FEARLESS run as compared to the standard run 
(Fig. [T2|) . This result is consistent both with the locally 
increased dissipation of turbulent to internal energy pro- 
vided by the SGS model and with the higher degree of 
mixing induced in the cluster core, shown by Urms in Ta- 
ble [1 

The radial distribution of turbulent energy is displayed 
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Fig. 14. — Same as Fig.|9] but showing the mass- weighted baryon 
radial velocity. 



TABLE 3 

Mass-weighted averages in a volume of 
[512 X 768 X 1280] kpc h"^, containing a subclump 

AND ITS wake, at 2 = 0.05. 



Quantity 


Adiabatic run 


fearless run 






1.13 


etot [10^*^ cm2 s-2] 


1.7447 


1.5746 


eint [lO^S cm2 s-2] 


5.4281 


5.9607 


fikin [10^® cm^ s-2] 


1.2032 


0.9786 


et [10" cm2 s-2] 




2.290 


Pt/(pt +Pthcrm) [%] 




3.70 



in Fig. [13] in terms of the turbulent velocity scaled to 
^min- The turbulent velocity at this length scale is below 
100 km s~^. There is a pronounced peak at r = 0.6 i?vir 
which is correlated with analogous trends in the turbu- 
lent pressure (Fig. [TT|) and in the radial velocity profile 
(Fig. [H)) . This structure is clearly linked with the most 
prominent merging clump shown in Fig. [8]i, ed analysed 
below in Table [3] 

There is an appealing similarity between the inter- 
vals of radii where q and the pressure ratio are larger 
(r < 0.1 i?vir and 0.4 i?vir < r < 0.8 i?vir), and the 
corresponding intervals where the temperature and en- 
tropy of the FEARLESS run are slightly larger than those 
computed for the adiabatic run. The opposite trend oc- 
curs in the interval in between, where Wturb is compar- 
atively smaller. The effects are very small, but suggest 
that the SGS model plays the same role in the ICM that 
was shown above for the cluster core, and in ! j5.3l for the 
global quantities. In case of radial profiles, the spheri- 
cal averaging combined with the intermittent behavior of 
turbulence tends to mask the turbulent effects. This can 
be better understood with a comparison of the values of 
q in Fig. [T2] and in the right-hand panels of Figs. [7] and 
[U the peak values in the slices are much larger than the 
spherical averages in the profile. 

The idea that locally the turbulence and its modeling 
can play a sizeable role is further corroborated by the 
data in Table [21 reporting the analysis at z = 0.05 of a 
small volume (512 x 768 x 1280 kpc h^-'^) that contain one 
of the clumps presented in i j5.2l and its wake (cf. Fig.[5]i). 
The morphology of this accreting subcluster in the fear- 



less run is not substantially different from the the adi- 
abatic one. The energy content, however, is rather dif- 
ferent from that in the cluster core: in the region under 
consideration ekin is dominant with respect to Cmt- The 
importance of the turbulence, injected by the hydrody- 
namical instabilities in the wake of the moving clump, 
is testified by the large ratio S/e and by the turbulent 
pressure support, which is at the level of some percent, 
about one order of magnitude larger than the spherical 
averages in Fig. [TlJ Despite of the slightly smaller av- 
erage of etot in the fearless run with respect to the 
adiabatic one, one can see an increase of eint, mostly at 
the expenses of ekin, resulting from the turbulent dissi- 
pation. The decrease of ekin is rather large (19%) but 
can be partially ascribed to the difhculty of comparing 
energy budgets in such open volumes. 

It is important to stress the deep difference between 
the turbulent velocity profile in Fig. [T3land the profile of 
Vrms, defined bv equation (j^ ) (see also lNorman fc BrvanI 
119991 : llapichino fc NiemeveH [2008 ^ . In the former case, 
the mass- weighted average of a local quantity (i.e., de- 
fined in every cell) is computed for each spherical shell, 
whereas in the latter interpreted shell-wise 

as the standard deviation with respect to the average 
(v). Clearly, the latter definition does not retain any 
information related to a length scale, and can be inter- 
preted as turbulent velocity only in a loose sense. From 
this point of view, the turbulent velocity provided by the 
SGS model is a more powerful probe of the features of a 
turbulent flow. On the other hand, spherically averaged 
velocity dispersions (and the derived turbulent pressure) 
are meaningful in comparison with observations, for ex- 
am ple in the procedu re for estimating the cluster mass 
(cf. iRasia et "all |2006[ ) . According to this different def- 
inition, the spherically averaged turbulent pressure of 
the simulated cluster (in a run similar to t he adiabatic 
one pr esented here) is reported in la pichino fc Niemeveii 
(|2008f ). It reaches values around 10%, in agree ment with 
the va lues found recently in simulations by iLau et ahl 
(|2009| ) for relaxed clusters. The turbulent pressure is 
somewhat smaller in the fearless run because of its 
slightly reduced content in kinetic energy, but the differ- 
ence is small, and is not expected to significantly affect 
the estimates of the cluster mass. 

6. DISCUSSION AND CONCLUSIONS 

Large Eddy Simulations (LES) are based on the no- 
tion of filtering the fiuid dynamic equations at a spe- 
cific length scale, thus performing a scale separation be- 
tween the resolved and the unresolved flow. The latter 
is treated by means of a subgrid scale model, which in 
turn is coupled to the hydrodynamical equations govern- 
ing the former. In principle, a single scale separation 
is incompatible with adaptive mesh refinement (AMR) 
codes, often used to study astrophysical phenomena. 

One of the aims of the present work was to address 
this numerical problem by means of developing, imple- 
menting, and applying a new numerical scheme that uses 
AMR and LES in combination, which we called fear- 
less. This novel tool is suitable for modeling turbulent 
flows over a wide range of length scales, a key feature in 
the treatment of many astrophysical flows including the 
intra-cluster medium. 

We showed that the idea of our approach to cor- 
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rect the velocity and kinetic energy at grid refine- 
ment /derefinement, according to local Kolmogorov scal- 
ing, produces consistent results in simulations of driven 
turbulence. We demonstrated that energy conservation 
and the scaling of turbulent energy in our adaptive sim- 
ulations is consistent with static grid simulations. 

To our knowledge, this work shows the first applica- 
tion of an SGS model to simulations of the formation and 
evolution of a galaxy cluster. The results give rise to sev- 
eral interesting implications with regard to the physics of 
galaxy clusters and to the numerical methods employed 
for their exploration in computational cosmology. 

The production of turbulence induced by minor 
merge rs, analytically studied by ISubramanian et al.l 
( 20061) and a ddressed by several numeri cal inve stigations 
jHeinz et al.|[2003: Takizaw a 2005a.bl: lAsai et al.. 2007;: 
iDursi fc Pfrommeii 2008: la pichino et al.ll2008D . is accu- 
rately tracked by the newly defined turbulent subgrid 
energy (Figs. [7] and [5]), although the level of resolution of 
the idealised setups cannot be reached by cosmological 
simulations. The visualization and subsequent analysis 
and postprocessing of turbulence and related quantities 
is therefore easier and more consistent. Turbulence in the 
ICM appears to be subsonic, in agreement with previous 
results. The average ratio between the dissipation and 
the production term E/e in the SGS model is close to 
unity, namely typical of a system where the turbulence 
is roughly stationary. On the other hand, this ratio is 
locally variable (cf. Table [2|) and, together with the in- 
terr nittent nature of tu rbulence in the ICM {^5^ see 
also llapichino &: Niemev cr 2008), delivers the picture of 
a fiow where turbulent motions are randomly ini tiated 
by merger events and the n gradually decay (.FrischlfTQOSl : 
ISubramanian et all 120061 ). We also notice that, in sim- 
ulations of driven turbulence in a periodic box (Fig. ^ , 
the decrease of ekin (noticed in our cluster simulation) is 
linked to an increase of et only in the early driving phase, 
not in the later equilibrium stage. 

The morphological evolution of the minor merger 
events and the subsequent injection of turbulence in the 
ICM (Figs. [7] and [8]) appear to be rather localized and 
intermittent, confirming the fe ature of turbulent flows as 
being not very volume-fillin g (jSubramanian et al.l 120061 : 
llapichino &: Niemeveill2008l ). The dissipation of turbu- 
lent to internal energy is thus modelled as a markedly lo- 
cal process, consistent with the theoretical expectations. 

The effect of the SGS model on the cluster energy bud- 
get is well exemplified by the comparison of our simula- 
tions at z = (Table [1]). Although the value of et is 
small compared to eint, this energy buffer is locally ef- 
fective in transferring the kinetic energy to the thermal 
component. The dissipative effects are therefore more 
relevant in those locations where et is relatively large, 
like the cluster core (Fig. [T3|) . In general, the main con- 
tribution of the FEARLESS approach is to add a more 
physically motivated contribution to the energy dissipa- 
tion, which in Eulerian codes is otherwise purely numer- 
ical. In FEARLESS, part of the energy flux from resolved 
scales to the thermal reservoir is retained in the buffer 
turbulent energy, et, and is further dissipated (turbulent 
dissipation) according to a local and temporal evolution 
determined by the SGS model. 

Besides local effects, the importance of the SGS model 
for the overall cluster structure appears small, because of 



the modest subgrid energy contribution (Fig. Ilip . One 
remark about the simulat ed cluster is important at thi s 
point: as also verifled in llapichino fc Niemevej (|2008l), 
this structure is very relaxed (see also Fig. [14)) . Simula- 
tions of more perturbed structures w ith recent or ongoing 
major mergers are in preparation ( Paul et al.|[2009 ). be- 
cause they will help to clarify the role of the turbulent 
energy (and of its modelling) in the cluster energy budget 
in cases where its magnitude is larger. From this view- 
point, the radial increase of turbulent pressure support 
in the cluster outskirt (Fig. [TT|) is interesting for physi- 
cal mechanisms (like the acceleration of cosmic rays and 
magnetic field amplification) where the knowledge of the 
turbulent state of the flow is needed. 

More turbulence in the cluster core is required, for ex- 
ample, to reproduce the iron abundance profi l e in c ool 
core clusters. Following iDennis fc ChandranI (|2005l ). a 
turbulent diffusion coefficient can be defined as Dturb ~ 
0.1 q I, where q is the turbulent velocity at the length 
scale I. Using / — laun and q 60 km s~^, we find 
£>turb — 2 X 10^^ cm^ s~^. We notice that this value is 
smaller than the estimates of the effective diffusion coef- 
ficient in the cluster models of iRebusco et al.l ()2005f ) and 
iRebusco et al.l ()2006D , which aim to reproduce the turbu- 
lent diffusion of metals in the cores of selected clusters. 
In particular, the cited models require much larger tur- 
bulent velocities. In the framework of our cluster simu- 
lation, these velocities could be injected into the ICM by 
a vigorous merger event. Another possibility, explicitly 
suggested by the authors cited above, is to invoke the 
action of an AGN outflow as an additional stirring agent 
in the cluster core. 

The enhanced temperature profile in the fearless 
run is somehow reminiscent of the theoretical predic- 
tions about the role of turbulent heating in cluster cores 
([Dennis & Chandran 200j). We notice an apparent mis- 
understanding in the literature regarding this point. In 
our model (and in the theory of turbulent fiows in gen- 
eral), the dissipation of turbulent energy does not act 
as an additional energy source but simply releases the 
energy arising from the virialization process on a longer 
timescale than the quick shock heating. Nevertheless, 
we showed that turbulence, and the turbulent dissipa- 
tion as well, can be rather localized. Naively, one could 
think that an effective turbulent heating in cool cores 
would require a peak of turbulent energy in the cluster 
core, whose existence and magnitude should be justified 
theoretically. Again, the stirring induced by AGN ac- 
tivity is an open possibility which deserves further in- 
vcstigation. However, the model of iDennis fc ChandranI 
(2005) includes radiative cooling and thermal conduc- 
tion, and a detailed comparison is beyond the scope of 
the present work. We observe that additional physics 
which is here not addressed (thermal conduction, mag- 
netic fields) could bring further interesting implications 
for t he energy budget in the ICM and the turbulent mix- 
ing (Shar ma et al.l l2009'. and references therein). 

Consequent to both the enhanced dissipation and fluid 
mixing is the larger value of entropy in the fearless 
cluster core. A long-standing problem in cluster simula- 
tions is the shape of the entropy profile, which smoothly 
decreases in the center in SPH simulations whereas it 
flattens inside the core in runs with grid-based codes 
(jFrenk et al.l I1999D . This issue has been debated re- 
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cently in several works (among others. [Polag et al1l2005l . 
IWadslev et all l2008l iKawata et all |2009D, because it is 
controversial which difference of SPH and mesh-based 
codes it results from. It has been claimed that the source 
of di screpancy prob ably lies in the treatment of fluid mix- 
ing ([Mitchell et al.ir2009) : the weaknesses of SPH in this 
regard are known, but the ability of mesh codes to model 
the turbulent cascade on length scales comparable with 
the grid resolution has not been addressed in a satisfac- 
tory way. It is therefore unclear whether the flat core 
entropy in grid codes correctly represents the physics of 
the ICM, or perhaps numerical effects harm th e robust- 
ness of this feature. Recently, ISpringell (|2009 t) pointed 
out that the core temperature and entropy in grid-based 
codes are affected by a spurious increase, caused by the 
N-body noise in the gravitational force field. In our opin- 
ion, the higher entropy core value in the fearless run 
suggests that the typical flat entropy core is a hydrody- 
namical feature which requires a better understanding of 
the numerics in mesh codes, and is at least not primarily 
caused by N-body noise. 

The SGS model applied in this work has to be con- 
sidered as an intermediate solution to address some 
basic questions related to dynamics of the turbulent 
intra-cluster medium. A more elaborate model that 



is able to handle the complexity of the flow (wide 
range of Mach numbers and large density gradients 
as well as pronounced inhomogeneities) in simulations 
of lar ge scale structure evoluti on is under develop- 
ment ( jSchmidt fc Federrathl[2009l ) . This first application 
shows the promising perspectives for the use of an SGS 
model in combination with AMR and its potential im- 
pact on many branches of numerical astrophysics. 
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